Обнаружение статистически значимых отличий в уровнях экспрессии генов больных раком

Это задание поможет вам лучше разобраться в методах множественной проверки гипотез и позволит применить ваши знания на данных из реального биологического исследования.

В этом задании вы:

вспомните, что такое t-критерий Стьюдента и для чего он применяется сможете применить технику множественной проверки гипотез и увидеть собственными глазами, как она работает на реальных данных почувствуете разницу в результатах применения различных методов поправки на множественную проверку

Основные библиотеки и используемые методы:

Библиотека scipy и основные статистические функции:http://docs.scipy.org/doc/scipy/reference/stats.html#statistical-functions

Библиотека statmodels для методов коррекции при множественном сравнении:

http://statsmodels.sourceforge.net/devel/stats.html

Статья, в которой рассматриваются примеры использования statsmodels для множественной проверки гипотез:

http://jpktd.blogspot.ru/2013/04/multiple-testing-p-value-corrections-in.html

Описание используемых данных

Данные для этой задачи взяты из исследования, проведенного в Stanford School of Medicine. В исследовании была предпринята попытка выявить набор генов, которые позволили бы более точно диагностировать возникновение рака груди на самых ранних стадиях.

В эксперименте принимали участие 24 человек, у которых не было рака груди (normal), 25 человек, у которых это заболевание было диагностировано на ранней стадии (early neoplasia), и 23 человека с сильно выраженными симптомами (cancer).

Ученые провели секвенирование биологического материала испытуемых, чтобы понять, какие из этих генов наиболее активны в клетках больных людей.

Секвенирование — это определение степени активности генов в анализируемом образце с помощью подсчёта количества соответствующей каждому гену РНК.

В данных для этого задания вы найдете именно эту количественную меру активности каждого из 15748 генов у каждого из 72 человек, принимавших участие в эксперименте.

Вам нужно будет определить те гены, активность которых у людей в разных стадиях заболевания отличается статистически значимо.

Кроме того, вам нужно будет оценить не только статистическую, но и практическую значимость этих результатов, которая часто используется в подобных исследованиях.

Диагноз человека содержится в столбце под названием "Diagnosis".

Практическая значимость изменения

Цель исследований — найти гены, средняя экспрессия которых отличается не только статистически значимо, но и достаточно сильно. В экспрессионных исследованиях для этого часто используется метрика, которая называется fold change (кратность изменения). Определяется она следующим образом:

$$F_{c}(C,T) = \begin{cases} \frac{T}{C}, T>C \\ -\frac{C}{T}, T<C \end{cases}$$

где C,T — средние значения экспрессии гена в control и treatment группах соответственно. По сути, fold change показывает, во сколько раз отличаются средние двух выборок.

Инструкции к решению задачи

Задание состоит из трёх частей. Если не сказано обратное, то уровень значимости нужно принять равным 0.05.

Часть 1: применение t-критерия Стьюдента

В первой части вам нужно будет применить критерий Стьюдента для проверки гипотезы о равенстве средних в двух независимых выборках. Применить критерий для каждого гена нужно будет дважды:

  • для групп normal (control) и early neoplasia (treatment)
  • для групп early neoplasia (control) и cancer (treatment)

В качестве ответа в этой части задания необходимо указать количество статистически значимых отличий, которые вы нашли с помощью t-критерия Стьюдента, то есть число генов, у которых p-value этого теста оказался меньше, чем уровень значимости.


In [31]:
import pandas as pd
import scipy.stats

In [32]:
df = pd.read_csv("gene_high_throughput_sequencing.csv")
control_df = df[df.Diagnosis == 'normal']
neoplasia_df = df[df.Diagnosis == 'early neoplasia']
cancer_df = df[df.Diagnosis == 'cancer']

In [33]:
# scipy.stats.ttest_ind(data.Placebo, data.Methylphenidate, equal_var = False)

In [34]:
genes = filter(lambda x: x not in ['Patient_id', 'Diagnosis'], df.columns.tolist())

In [35]:
control_vs_neoplasia = {}
neoplasia_vs_cancer = {}

In [36]:
for gene in genes:
    control_vs_neoplasia[gene] = scipy.stats.ttest_ind(control_df[gene], neoplasia_df[gene], equal_var = False).pvalue
    neoplasia_vs_cancer[gene] = scipy.stats.ttest_ind(cancer_df[gene], neoplasia_df[gene], equal_var = False).pvalue

In [37]:
print control_df['LOC643837'],neoplasia_df['LOC643837']


0     1.257614
1     4.567931
2     2.077597
3     2.066576
4     2.613616
5     3.942275
6     1.084113
7     3.153900
8     2.551800
9     3.693128
10    3.558222
11    0.938061
12    1.003451
13    7.364879
14    2.561871
15    2.971205
16    2.871767
17    1.045382
18    1.801865
19    3.515834
20    2.234576
21    4.717822
22    1.474174
23    1.282995
Name: LOC643837, dtype: float64 24    2.516305
25    1.937270
26    1.405858
27    2.131757
28    2.421766
29    4.668232
30    3.386331
31    1.247440
32    1.591747
33    0.979074
34    4.730638
35    0.964837
36    0.949252
37    1.068210
38    1.679849
39    0.885375
40    4.567425
41    3.060481
42    2.888731
43    1.285190
44    5.753165
45    3.983593
46    4.391919
47    1.122588
48    3.155315
Name: LOC643837, dtype: float64

In [38]:
scipy.stats.ttest_ind(control_df['LOC643837'], neoplasia_df['LOC643837'], equal_var = False).pvalue


Out[38]:
0.69076601574973551

In [39]:
control_vs_neoplasia_df = pd.DataFrame.from_dict(control_vs_neoplasia, orient = 'index')
control_vs_neoplasia_df.columns = ['control_vs_neoplasia_pvalue']

In [40]:
neoplasia_vs_cancer_df = pd.DataFrame.from_dict(neoplasia_vs_cancer, orient = 'index')
neoplasia_vs_cancer_df.columns = ['neoplasia_vs_cancer_pvalue']

In [42]:
neoplasia_vs_cancer_df


Out[42]:
neoplasia_vs_cancer_pvalue
UBE2Q1 0.453480
HIF3A 0.000654
UBE2Q2 0.114018
RNF10 0.702235
RNF11 0.346222
LRRC31 0.419521
RNF13 0.485952
REM1 0.003971
REM2 0.698901
C16orf13 0.049652
MZT2A 0.072839
MZT2B 0.163551
NDP 0.872778
PMM2 0.014564
PMM1 0.562594
ASS1 0.767728
NCBP1 0.886632
ZNF709 0.000414
ZNF708 0.151118
ZNF879 0.651241
CBX3 0.024891
NID2 0.501976
CAMK1 0.729918
STYXL1 0.059709
ZNF701 0.874405
ZNF700 0.151003
ZNF707 0.023744
LOC147670 0.820141
ZNF704 0.092377
ZC3H10 0.912989
... ...
PTRF 0.043195
LUZP6 0.942831
C1orf198 0.194484
CFL2 0.584340
CFL1 0.068066
C8orf83 0.669453
C1orf192 0.014639
PLAA 0.498337
C1orf190 0.004657
BCL6B 0.009442
MRPL23 0.279824
C1orf194 0.201948
SELL 0.038052
PLEKHG3 0.091142
NFIA 0.074839
SELO 0.131847
PLEKHG6 0.438145
PLEKHG4 0.002692
PLEKHG5 0.128011
SELE 0.488034
SLC7A10 0.024658
C2orf71 0.071053
SIGLEC1 0.000088
AP4M1 0.624477
GNGT2 0.003862
AIP 0.909838
SERPINH1 0.006149
NFIX 0.143400
SELP 0.812087
SELS 0.131701

15748 rows × 1 columns


In [43]:
pvalue_df = control_vs_neoplasia_df.join(neoplasia_vs_cancer_df)
pvalue_df.head()


Out[43]:
control_vs_neoplasia_pvalue neoplasia_vs_cancer_pvalue
UBE2Q1 0.908853 0.453480
HIF3A 0.115963 0.000654
UBE2Q2 0.355621 0.114018
RNF10 0.081252 0.702235
RNF11 0.438372 0.346222

In [44]:
pvalue_df[pvalue_df.control_vs_neoplasia_pvalue < 0.05].shape


Out[44]:
(1575, 2)

In [45]:
pvalue_df[pvalue_df.neoplasia_vs_cancer_pvalue < 0.05].shape


Out[45]:
(3490, 2)

Часть 2: поправка методом Холма

Для этой части задания вам понадобится модуль multitest из statsmodels.

import statsmodels.stats.multitest as smm

В этой части задания нужно будет применить поправку Холма для получившихся двух наборов достигаемых уровней значимости из предыдущей части. Обратите внимание, что поскольку вы будете делать поправку для каждого из двух наборов p-value отдельно, то проблема, связанная с множественной проверкой останется.

Для того, чтобы ее устранить, достаточно воспользоваться поправкой Бонферрони, то есть использовать уровень значимости 0.05 / 2 вместо 0.05 для дальнейшего уточнения значений p-value c помощью метода Холма.

В качестве ответа к этому заданию требуется ввести количество значимых отличий в каждой группе после того, как произведена коррекция Холма-Бонферрони. Причем это число нужно ввести с учетом практической значимости: посчитайте для каждого значимого изменения fold change и выпишите в ответ число таких значимых изменений, абсолютное значение fold change которых больше, чем 1.5.

Обратите внимание, что

применять поправку на множественную проверку нужно ко всем значениям достигаемых уровней значимости, а не только для тех, которые меньше значения уровня доверия. при использовании поправки на уровне значимости 0.025 меняются значения достигаемого уровня значимости, но не меняется значение уровня доверия (то есть для отбора значимых изменений скорректированные значения уровня значимости нужно сравнивать с порогом 0.025, а не 0.05)!


In [46]:
import statsmodels.stats.multitest as smm

In [47]:
pvalue_df['control_mean_expression'] = control_df.mean()
pvalue_df['neoplasia_mean_expression'] = neoplasia_df.mean()
pvalue_df['cancer_mean_expression'] = cancer_df.mean()

In [48]:
def abs_fold_change(c, t):
    if t > c:
        return t/c
    else:
        return c/t

In [49]:
pvalue_df['control_vs_neoplasia_fold_change'] = map(lambda x, y: abs_fold_change(x, y), 
                                                    pvalue_df.control_mean_expression,
                                                    pvalue_df.neoplasia_mean_expression
                                                   )
pvalue_df['neoplasia_vs_cancer_fold_change'] = map(lambda x, y: abs_fold_change(x, y),
                                                   pvalue_df.neoplasia_mean_expression,
                                                   pvalue_df.cancer_mean_expression
                                                  )

In [50]:
pvalue_df['control_vs_neoplasia_rej_hb'] = smm.multipletests(pvalue_df.control_vs_neoplasia_pvalue, alpha=0.025, method='h')[0]
pvalue_df['neoplasia_vs_cancer_rej_hb'] = smm.multipletests(pvalue_df.neoplasia_vs_cancer_pvalue, alpha=0.025, method='h')[0]

In [51]:
pvalue_df[(pvalue_df.control_vs_neoplasia_rej_hb) & (pvalue_df.control_vs_neoplasia_fold_change > 1.5)].shape


Out[51]:
(2, 9)

In [52]:
pvalue_df[(pvalue_df.neoplasia_vs_cancer_rej_hb) & (pvalue_df.neoplasia_vs_cancer_fold_change > 1.5)].shape


Out[52]:
(77, 9)

Часть 3: поправка методом Бенджамини-Хохберга

Данная часть задания аналогична второй части за исключением того, что нужно будет использовать метод Бенджамини-Хохберга.

Обратите внимание, что методы коррекции, которые контролируют FDR, допускает больше ошибок первого рода и имеют большую мощность, чем методы, контролирующие FWER. Большая мощность означает, что эти методы будут совершать меньше ошибок второго рода (то есть будут лучше улавливать отклонения от H0, когда они есть, и будут чаще отклонять H0, когда отличий нет).

В качестве ответа к этому заданию требуется ввести количество значимых отличий в каждой группе после того, как произведена коррекция Бенджамини-Хохберга, причем так же, как и во второй части, считать только такие отличия, у которых abs(fold change) > 1.5.


In [57]:
pvalue_df['control_vs_neoplasia_rej_bh'] = smm.multipletests(pvalue_df.control_vs_neoplasia_pvalue, alpha=0.025, method='fdr_i')[0]
pvalue_df['neoplasia_vs_cancer_rej_bh'] = smm.multipletests(pvalue_df.neoplasia_vs_cancer_pvalue, alpha=0.025, method='fdr_i')[0]

In [60]:
pvalue_df.control_vs_neoplasia_rej_bh.value_counts()


Out[60]:
False    15744
True         4
Name: control_vs_neoplasia_rej_bh, dtype: int64

In [54]:
pvalue_df[(pvalue_df.control_vs_neoplasia_rej_bh) & (pvalue_df.control_vs_neoplasia_fold_change > 1.5)].shape


Out[54]:
(4, 11)

In [55]:
pvalue_df[(pvalue_df.neoplasia_vs_cancer_rej_bh) & (pvalue_df.neoplasia_vs_cancer_fold_change > 1.5)].shape


Out[55]:
(524, 11)

In [ ]: